/*
    Linear Programming Library
    
    Created by Team Paragon (Shena Hoffmann, Arthur Powalka, Nasir Kenarangui)
    CSCE 689
    
    Karmarkar.cpp
    Defines the variables and functions used for providing functionality related to the Karmarkar linear programming algorithm
*/

#include "Karmarkar.h"

// ............ CONSTRUCTORS ............ //

Karmarkar::Karmarkar()
{ //empty constructor 
}

Karmarkar::Karmarkar(const LinProg& lp):LinProg(lp) //dispatches lp to LinProg copy constructor
{ }

Karmarkar::Karmarkar(const Karmarkar& lp):LinProg(lp) //dispatches lp to LinProg copy constructor
{ }

Karmarkar::Karmarkar(bool max, vector<double> objective, vector<vector<double>> constraints, vector<string> inequalitySigns, vector<double> inequalityVals, vector<bool> constrainedVars)
{
    this->max = max;
    this->objective = objective;
    this->constraints = constraints;
    this->inequalitySigns = inequalitySigns;
    this->inequalityVals = inequalityVals;
    this->constrainedVars = constrainedVars;
    verifyLinProg();
}

// ............ DESTRUCTOR ............ //
Karmarkar::~Karmarkar()
{

}

// ............ HELPER FUNCTIONS ............ //

mat Karmarkar::improvement(mat &x_k, mat &c, mat &A, mat &e)
{		
	double n=e.n_rows;
	double alpha=(n-1)/(3*n);
	double r=1/sqrt(n*(n-1));
	mat X_k=diagmat(x_k);	//X_k and x_k are different
	mat B=join_cols(A*X_k,trans(e));		
	mat temp=trans(B)*inv(B*trans(B))*B;
	mat temp2(temp.n_rows,temp.n_cols);
	temp2.eye();
	mat P_B=temp2-temp;
	mat delta_x=(-1)*P_B*X_k*c;
	mat y_k1=(e/n)+alpha*(1/sqrt(n*(n-1)))*delta_x/norm(delta_x,2);
	x_k=(X_k*y_k1)/as_scalar(trans(e)*X_k*y_k1);
	return x_k;
}

double Karmarkar::findUpperBound()
{
	double upperBound=0;
	double temp;
	for(int h=0;h<constraintsA.n_rows;h++)	
	{
		for_each(constraintsA.begin_row(h),constraintsA.end_row(h),[&](double denom)
		{
			if (denom!=0)	//There was a huge logic error here. Previously it was denom>0.... big mistake
			{
				temp=valuesB[h]/denom;
				
				if(temp > upperBound)
					upperBound=temp;
			}	
		});
	}
	return upperBound;
}

bool Karmarkar::isFeasible(const mat &solutionMatrix)
{
	bool feasible=true; //initially considered feasible
	colvec lhs=constraintsA*solutionMatrix; //evaluates the constraints at the solution found to check if they are in the feasible region
	double hundred=100;
	for(int t=0;t<lhs.n_elem;t++) //compares the solution to the constraintValues
	{
		if( !(lhs(t)<=valuesB(t)) )
		{
			//if( !(abs(lhs(t)-valuesB(t))<.1) )	//make tolerance higher. needs to be a percentage of solution	
			if( !(abs(lhs(t)-valuesB(t))/hundred <0.05 )) //needs to be tolerant because solution may be slightly smaller/larger than optimal		
				feasible=false;		        //since they are not exact they may break feasibility by a tiny margin
		}						
	}	
	return feasible;
}

void Karmarkar::setClassVariables(LinProg &obj) //if we make argument Karmarkar then we will have direct access to protected variables
{
	rowvec temp = conv_to<rowvec>::from(obj.getObjective());
	objectiveC=temp;
	vector< vector<double> > conVec=obj.getConstraints();	//obj's type doesn't know about being Karmarkar's base class, so we don't have direct access to protected variables.

	for_each(conVec.begin(),conVec.end(), [&](vector<double> value)
	{
		rowvec temp = conv_to<rowvec>::from(value);
		constraintsA=join_cols(constraintsA,temp);
	});

	valuesB = conv_to<mat>::from(obj.getInequalityVals());	
}

// ............ FUNCTIONS FOR THE USER ............ //

Karmarkar Karmarkar::karForm(mat C, mat A, mat b)
{	
//cout<<C<<endl;
//cout<<A<<endl;
//cout<<b<<endl;
//cout<<endl;	
	verifyLinProg();
	mat cDual=trans(b);
	mat aDual=trans(A);	//constraint coefficients for the dual
	mat bDual=trans(C);	//b values for the dual

//cout<<cDual<<endl;
//cout<<aDual<<endl;
//cout<<bDual<<endl;
//cout<<endl;

	double rows=aDual.n_rows;
	double cols=aDual.n_cols;
	double r=A.n_rows;
	double c=A.n_cols;
	mat slackP(r,r);
	slackP.eye();
	mat slackD(rows,rows);
	slackD.eye();
	slackD=(-1)*slackD;	
	A=join_rows(A,slackP);
	aDual=join_rows(aDual,slackD);
	
	rows=aDual.n_rows;
	cols=aDual.n_cols;
	r=A.n_rows;
	c=A.n_cols;
	
	A=join_rows(A,zeros(r,cols));
	aDual=join_rows(zeros(rows,c),aDual);
	
	rows=aDual.n_rows;
	cols=aDual.n_cols;
	r=A.n_rows;
	c=A.n_cols;
	
	mat bCon=ones(1,cols+1);
	A=join_rows(A,zeros(r,1));
	aDual=join_rows(aDual,zeros(rows,1));
	varBound=findUpperBound()*cols;	//findUpperBound will find an upperbound on the solution. This is required in order to put linear problem into Karmarkar form
	mat bbCon={varBound};
	A=join_rows(A,-1*b);
	
	aDual=join_rows(aDual,-1*bDual);
	bCon=join_rows(bCon,-1*bbCon);
	double k=slackP.n_cols;
	double m=bCon.n_cols;
	
	mat aFinal=join_rows(C,zeros(1,k));
	aFinal=join_rows(aFinal,-1*cDual);

	double v=aFinal.n_cols;	
	aFinal=join_rows(aFinal,zeros(1,m-v));
	aFinal=join_cols(aFinal,A);
	aFinal=join_cols(aFinal,aDual);
	aFinal=join_cols(aFinal,bCon);

	r=aFinal.n_rows;

	mat extra_col=sum(aFinal,1);	//1 indicates sum of each row
	aFinal=join_rows(aFinal,-1*extra_col);
	c=aFinal.n_cols;
	mat obj=zeros(1,c);
	obj(c-1)=1;
	obj=trans(obj);
	mat e=ones(1,c);
	aFinal=join_cols(aFinal,e);
	Karmarkar karStruct;	
	karStruct.objectiveC=obj;
	karStruct.constraintsA=aFinal;
//cout<<endl<<aFinal<<endl;
	return  karStruct;
}

tuple<vector<double>, string, double> Karmarkar::solve()
{
	LinProg standardForm=convertToStandardForm(); //return the standard form of this* object
	//standardForm.printLinProg();
	setClassVariables(standardForm); //puts the standard form of problem into mat objectiveC, mat constraintsA, mat valuesB
//cout<<objectiveC<<endl;
//cout<<valuesB<<endl;
//cout<<constraintsA<<endl;
	Karmarkar first=karForm(objectiveC,constraintsA,valuesB); //creates the karmarkar form of the problem
//cout<<first.constraintsA<<endl;
	int f=first.constraintsA.n_rows-1;
	mat e=first.constraintsA.row(f);
//cout<<e<<endl;
	first.constraintsA.shed_row(f);

//cout<<first.constraintsA<<endl;
	e=trans(e);		
	double n=e.n_rows;
	mat x_0=e/n; 
	mat temp=improvement(x_0,first.objectiveC,first.constraintsA,e);
//cout<<(varBound+1)*temp<<endl;
//cout<<first.objectiveC<<endl;
	while(abs(as_scalar(trans(first.objectiveC)*temp))>.00001)
	{
		temp=improvement(temp,first.objectiveC,first.constraintsA,e);
//cout<<(varBound+1)*temp<<endl;
	}
	
//cout<<(varBound+1)*temp<<endl;
	/////////////FORMING SOLUTION TUPLE BELOW//////////////////////		
	//temp.set_size(objectiveC.n_cols,1);
	temp.reshape(objectiveC.n_cols,1,1);
//cout<<(varBound+1)*temp<<endl;
	mat truncatedSolution=(varBound+1)*temp;
	vector<double> varSolution;
	for(int f=0;f<truncatedSolution.n_rows;f++) //creates solution vector for the variables
	{
		varSolution.push_back(truncatedSolution(f));
	}

	mat zValue=objectiveC*truncatedSolution;
	double objValue=zValue(0,0);
	
	if(max==0)
	{
		objValue=(-1)*objValue;
	}
	string feasible;
	if (isFeasible(truncatedSolution))
	{
		feasible="feasible";
	}
	else
	{
		feasible="not feasible";
		varSolution.clear();
		objValue=0;
	}
	//Karmarkar does not solve unbounded problems... this is a possible check for unboundedness
	bool unBounded=0;	
	if(feasible=="feasible")	
	
	{
		for(int i=0;i<objectiveC.n_cols;i++)
		{
			if(truncatedSolution[i]<.00001)
			{
				//truncatedSolution[i]=0;
				unBounded=1;
			}
			else
				unBounded=0;
		}
	}	
	
	if(unBounded==1)
	{
		feasible="possibly unbounded";
		varSolution.clear();
		objValue=0;
	}
	objectiveC.clear();
	constraintsA.clear();
	valuesB.clear();
	///////////////////////////////////////////////////////////////////////////////////////////

	auto solutionData=make_tuple(varSolution, feasible, objValue);
	
	return solutionData;
}

mat Karmarkar::getObjectiveC()
{
	return objectiveC; 
}		

mat Karmarkar::getConstraintsA()
{
	return constraintsA; 
}
	
mat Karmarkar::getValuesB()
{
	return valuesB;
}

// ............ OVERLOADED OPERATORS ............ //

Karmarkar Karmarkar::operator=(const Karmarkar &lp)
{
    max = lp.max;
    objective = lp.objective;
    constraints = lp.constraints;
    inequalitySigns = lp.inequalitySigns;
    inequalityVals = lp.inequalityVals;
    constrainedVars = lp.constrainedVars;
}


